---
title: 'TSMOs and Democratic Diffusion Analysis: Appendix With Full Analysis'
author: "Author Blinded for Peer Review"
output:
  pdf_document:
    toc: true
    toc_depth: 3
    latex_engine: xelatex
header-includes:
- \usepackage{setspace}
- \usepackage{rotating}
fontsize: 10pt
---

## Introduction

This appendix contains the full results for all robustness checks we conducted for the paper: "Activists Against Autocrats: TSMOs and Democratic Diffusion." Note that, while in the main text of the article we only report regression results for our primary dependent variables of interest (Polyarchy, Participatory Democracy, Free and Fair Elections, and Freedom of Association) in all tables in the Appendix we report results for all five of the Variety of Democracy Project's "top-level" democracy indexes and the five components of the Polyarchy score. This gives greater insight into the diffusing effects of TSMOs and IGOs. In most models, our measures of TSMO diffusion are not only significant predictors of participatory democracy and freedom of association, but also significant predictors of deliberative democracy and free expression. We consider these results to be further support for our hypothesis about TSMO diffusion being focused on more participatory and expressive elements of democracy, and less focused on institutional elements of democracy. 

```{r setup, results='hide', message=FALSE, echo=FALSE}

## Setup Chunk - Load necessary packages, workspace, and functions.

knitr::opts_chunk$set(echo = FALSE)
packs <- c("tidyverse","multiwayvcov", "stargazer","effects","car", 
           "magrittr","broom","dotwhisker","lmtest","texreg","sandwich") # Load Necessary Packages
lapply(packs,library,character.only = T) 

source("functions.R")

# Load Datasets created in data_prep script

all.working.datasets <- readRDS("all_working_datasets.rds") 

working.data.all.tsmos <- all.working.datasets[[1]]

working.data.hr.only <- all.working.datasets[[2]]

working.data.women.only <- all.working.datasets[[3]]

rm(all.working.datasets)

```

```{r wtd.lead, results='hide',warning=FALSE,message=FALSE}

attach(working.data.all.tsmos)

wtd.lead.mods <- tsmomodel(dvform = "lead",
                           ivform = "wtd",
                           controlvars = c("ioscore","dem.perf","cw","log_gdppc","neighbor","log.tsmo.mems"))

clustered <- map(wtd.lead.mods, ~ cluster.ses(.,clust = ccode))

detach(working.data.all.tsmos)

clustered <- map(clustered, function(x) {rownames(x)[2:5] <- c("tsmoscore","ioscore","dem.econ.perf","neighbors")
x })

```



```{r wtd.lead (top5), results='asis',warning=FALSE}

stargazer(clustered[c(1:5)],
          omit = "ccode",
          column.sep.width = "5pt",
          column.labels = c("Polyarchy","Liberal Dem.","Participatory Dem.","Deliberative Dem.","Egalitarian Dem"),
          covariate.labels = c("TSMO Diffusion","IO Score","Democracy Econ. Perf.","Geog. Diffusion","Cold War","GDP per cap (log","Number of TSMOs (log)"),
          star.cutoffs = c(0.05,0.01,0.001),
          title = "TSMO Diffusion Across Dimensions of Democracy",
          header = F
)

```

```{r wtd.lead (components), results='asis', warning=FALSE}

stargazer(clustered[c(6:10)],
          omit = "ccode",
          column.labels = c("Elected Officials","Free Elections","Free Association","Free Expression","Suffrage"),
          covariate.labels = c("TSMO Diffusion","IO Score","Democracy Econ. Perf.","Geog. Diffusion","Cold War","GDP per cap (log","Number of TSMOs (log)"),
          star.cutoffs = c(0.05,0.01,0.001),
          notes = "Fixed Effects Omitted, Robust St. Err. in Parentheses, * p < 0.05, ** p < 0.01,*** p < 0.001",
          notes.append = FALSE,
          notes.align = "l",
          title = "TSMO Diffusion Across Components of Polyarchy",
          header = F
)
```

\clearpage

## Different Model Specifications

The models below replicate our primary models while varying the form of the independent variable, dependent variable, and both. Our alternate independent variable follows replicates Pevehouse's procedure for measuring democratic diffusion through IGOs. That is to say, for each country $i$ we average the level of democracy in each of their TSMO networks and then assign the country the maximum out of these averages. Our alternate dependent variable is democracy at $t$ minus democracy at $t - 1$ (the first difference in democracy scores).

### Models with Democracy Level Dependent Variable and TSMOScore IV

```{r lead.score.new, results='asis',warning = FALSE}
attach(working.data.all.tsmos)
score.lead.mods <- tsmomodel(dvform = "lead",
                           ivform = "score",
                           controlvars = c("ioscore","dem.perf","cw","log_gdppc","neighbor","log.tsmo.mems"))

clustered <- map(score.lead.mods, ~ cluster.ses(.,clust = ccode))

detach(working.data.all.tsmos)

clustered <- map(clustered, function(x) {rownames(x)[2:5] <- c("tsmoscore","ioscore","dem.econ.perf","neighbors")
x })


stargazer(clustered[c(1,3,7:9)],
          star.cutoffs = c(0.05,0.01,0.001),
          column.labels = c("Polyarchy","Participatory Dem.","Free Elections","Free Association","Free Expression"),
          omit = "ccode",
          keep.stat = c("n","rsq"),
          covariate.labels = c("TSMO Diffusion (score)","IO Score","Democracy Econ. Perf.","Geog. Diffusion","Cold War","GDP per cap (log)","TSMO Memberships (log)"),
          header = FALSE,
          column.sep.width = "0pt",
          title = "TSMOScore Diffusion IV - Democracy Level DV",
          notes = "Fixed Effects Omitted, Robust St. Err. in Parentheses, * p < 0.05, ** p < 0.01,*** p < 0.001",
          notes.append = FALSE,
          notes.align = "l")


```

\clearpage

### Models with First Difference Dependent Variable 

```{r diff.wtd.new, results='asis', warning=FALSE}

attach(working.data.all.tsmos)
wtd.diff.mods <- tsmomodel(dvform = "diff",
                                 ivform = "wtd",
                                 controlvars = c("ioscore","dem.perf","cw","log_gdppc","neighbor","log.tsmo.mems"))

clustered <- map(wtd.diff.mods, ~ cluster.ses(.,clust = ccode))

detach(working.data.all.tsmos)

clustered <- map(clustered, function(x) {rownames(x)[2:5] <- c("tsmoscore","ioscore","dem.econ.perf","neighbors")
x })


stargazer(clustered[c(1,3,7:9)],
          star.cutoffs = c(0.05,0.01,0.001),
          column.labels = c("$\\Delta$ Polyarchy","$\\Delta$ Participatory","$\\Delta$ Free Elections","$\\Delta$ Free Association","$\\Delta$ Free Expression"),
          omit = "ccode",
          keep.stat = c("n","rsq"),
          covariate.labels = c("TSMO Diffusion","IO Score","Democracy Econ. Perf.","Geog. Diffusion","Cold War","GDP per cap (log)","TSMO Memberships (log)"),
          header = FALSE,
          column.sep.width = "0pt",
          title = "Spatial Lag TSMO Diffusion IV - First Difference DV",
          notes = "Fixed Effects Omitted, Robust St. Err. in Parentheses, * p < 0.05, ** p < 0.01,*** p < 0.001",
          notes.append = FALSE,
          notes.align = "l")

```

\clearpage

### Models with TSMOScore Independent Variable and First Difference Dependent Variable

```{r score.diff.new, results='asis',warning = FALSE}
attach(working.data.all.tsmos)
score.diff.mods <- tsmomodel(dvform = "diff",
                           ivform = "score",
                           controlvars = c("ioscore","dem.perf","cw","log_gdppc","neighbor","log.tsmo.mems"))

clustered <- map(score.diff.mods, ~ cluster.ses(.,clust = ccode))

detach(working.data.all.tsmos)

clustered <- map(clustered, function(x) {rownames(x)[2:5] <- c("tsmoscore","ioscore","dem.econ.perf","neighbors")
x })


stargazer(clustered[c(1,3,7:9)],
          star.cutoffs = c(0.05,0.01,0.001),
          column.labels = c("$\\Delta$ Polyarchy","$\\Delta$ Participatory","$\\Delta$ Free Elections","$\\Delta$ Free Association","$\\Delta$ Free Expression"),
          omit = "ccode",
          keep.stat = c("n","rsq"),
          covariate.labels = c("TSMO Diffusion (score)","IO Score","Democracy Econ. Perf.","Geog. Diffusion","Cold War","GDP per cap (log)","TSMO Memberships (log)"),
          header = FALSE,
          column.sep.width = "0pt",
          title = "TSMOScore Diffusion IV - Democracy Level DV",
          notes = "Fixed Effects Omitted, Robust St. Err. in Parentheses, * p < 0.05, ** p < 0.01,*** p < 0.001",
          notes.append = FALSE,
          notes.align = "l")

```

\clearpage

### Models Run Only on Non-Democracies

The following models replicate the modeling strategy that we report in the main text, but remove all country years with a polyarchy score above 0.5 (a rough threshold value for democracy, per the V-Dem Codebook). The results for TSMO diffusion are much stronger for Polyarchy and its sub-dimensions, while the measure of TSMO diffusion of participatory democracy loses significance (though it remains signed correctly). 

```{r nondems, results='asis',warning=FALSE}
nondems <- working.data.all.tsmos %>% filter(polyarchy < 0.5)

attach(nondems)

nondem.mods <- tsmomodel(dvform = "lead",
                         ivform = "wtd",
                         controlvars = c("ioscore","dem.perf","cw","log_gdppc","neighbor","log.tsmo.mems"))

clustered <- map(nondem.mods, ~ cluster.ses(.,clust = ccode))

detach(nondems)

clustered <- map(clustered, function(x) {rownames(x)[2:5] <- c("tsmoscore","ioscore","dem.econ.perf","neighbors")
x })


stargazer(clustered[c(1,3,7:9)],
          star.cutoffs = c(0.05,0.01,0.001),
          column.labels = c("Polyarchy","Participatory Dem.","Free Elections","Free Association","Free Expression"),
          omit = "ccode",
          keep.stat = c("n","rsq"),
          covariate.labels = c("TSMO Diffusion","IO Score","Democracy Econ. Perf.","Geog. Diffusion","Cold War","GDP per cap (log)","TSMO Memberships (log)"),
          header = FALSE,
          column.sep.width = "0pt",
          title = "Model Run on Non-Democratic Sub-Sample",
          notes = "Fixed Effects Omitted, Robust St. Err. in Parentheses, * p < 0.05, ** p < 0.01,*** p < 0.001",
          notes.append = FALSE,
          notes.align = "l")
```

\clearpage

### Main Models with Additional Controls

These models replicate the independent and dependent variable set-up in the primary models, but add several additional control variables. First, we include controls for time trends. These are simply the number of years since the beginning of the observation period (1953) plus that number's square and cube. We also include a measure of logged population from the UN, a measure of per capita oil and natural gas rents from Ross 2015, a logged measure of a country's total number of TSMO memberships, relative to the global median, and the global average polyarchy score as a rough measure of the overall level of democracy globally.

The main results from our primary models remain robust to the inclusion of these additional controls.

```{r controls, results='asis', warning = FALSE}
### Add Time Trend Transformation

more.controls <- working.data.all.tsmos %>%
  mutate(time.trend = year - 1953,
         time.trend.sq = time.trend^2,
         time.trend.cub = time.trend^3) %>% 
  group_by(year) %>% 
  mutate(avg_polyarchy = mean(polyarchy,na.rm = T)) %>% 
  ungroup()

attach(more.controls)

control.mods <- tsmomodel(dvform = "lead",
                           ivform = "wtd",
                           controlvars = c("ioscore","dem.perf","cw","log_gdppc","neighbor","log_pop","log_fuel","time.trend","time.trend.sq","time.trend.cub","log.tsmo.mems","avg_polyarchy"))

clustered <- map(control.mods, ~ cluster.ses(.,clust = ccode))

detach(more.controls)

clustered <- map(clustered, function(x) {rownames(x)[2:5] <- c("tsmoscore","ioscore","dem.econ.perf","neighbors")
x })

stargazer(clustered[c(1,3,7:9)],
          star.cutoffs = c(0.05,0.01,0.001),
          column.labels = c("Polyarchy","Participatory","Free Elections","Free Association","Free Expression"),
          keep.stat = c("n","rsq"),
          covariate.labels = c("TSMO Diffusion","IO Score","Democracy Econ. Perf.","Geog. Diffusion","Cold War","GDP per cap (log)", "Population (log)", "Fuel Rents per cap (log)","Time","Time (sq)","Time (cub)","TSMO Memberships (log)","Average Polyarchy"),
          header = FALSE,
          omit = "ccode",
          column.sep.width = "0pt",
          title = "Main Models With Extra Controls",
          notes = "Fixed Effects Omitted, Robust St. Err. in Parentheses, * p < 0.05, ** p < 0.01,*** p < 0.001",
          notes.append = FALSE,
          notes.align = "l")
```

\clearpage

### Models With Longer TSMO Lags

One concern in any observational study is reverse causality. In our case, the concern is that periods of nascent democratic change may attract more TSMOs. One way of addressing this concern is to see whether high TSMO scores have a democratic diffusion effect over the long-term, and do not just affect short-term democratic change. Thus we replicated our primary models, but extended the lag period between the point at which we measured the TSMO score and the change in democracy. We re-ran our primary models for lag periods *t* - 1 up to *t* - 10.  The figure below reports the results of these tests for the five independent variables that we focus on in the main text: the TSMO scores for Polyarchy, Participatory Democracy, Free and Fair Elections, Free Expression, and Freedom of Association. Our three hypotheses would suggest that TSMO diffusion variables for Polyarchy and Free and Fair Elections should have weak, insignificant effect, while the TSMO diffusion variables for Participatory Democracy, Freedom of Expression, and Freedom of Association should have larger, more robust effects. The points in the figure are the point estimate for the relevant TSMOscore measure, while the bars are a 95% confidence interval.

The most robust of our findings is the relationship with Freedom of Expression, which remains large and highly significant across all ten lags. The Participatory Democracy, Freedom of Association, and Free and Fair Elections TSMO diffusion variables retain significance up until *t* - 3, while the TSMO Polyarchy diffusion variable loses significance at *t* - 2.


```{r longer_lags, results='asis', warning=FALSE}

attach(working.data.all.tsmos)
lag1.mods <- tsmomodel(dvform = "lead",
                       ivform = "wtd",
                       controlvars = c("ioscore","dem.perf","cw","log_gdppc","neighbor","log.tsmo.mems"))
 ses <- lapply(lag1.mods, function(x) {
    sqrt(diag(cluster.vcov(x, ccode)))
  })
detach(working.data.all.tsmos)

lag1 <- bind_rows(tidy(lag1.mods[[1]])[2,],
                  tidy(lag1.mods[[3]])[2,],
                  tidy(lag1.mods[[7]])[2,],
                  tidy(lag1.mods[[8]])[2,],
                  tidy(lag1.mods[[9]])[2,]) %>%
  mutate(lag = 1,
         std.error = c(ses[[1]][2],
                       ses[[3]][2],
                       ses[[7]][2],
                       ses[[8]][2],
                       ses[[9]][2]))

long.lags.2 <- working.data.all.tsmos %>%
  group_by(ccode) %>%
  arrange(year) %>%
  mutate_at(vars(tsmo.wtd.polyarchy:tsmo.wtd.suffr), ~ lag(.,n = 1L)) %>%
  ungroup()

attach(long.lags.2)
lag2.mods <- tsmomodel(dvform = "lead",
                       ivform = "wtd",
                       controlvars = c("ioscore","dem.perf","cw","log_gdppc","neighbor","log.tsmo.mems"))
 ses <- lapply(lag2.mods, function(x) {
    sqrt(diag(cluster.vcov(x, ccode)))
  })
detach(long.lags.2)

lag2 <- bind_rows(tidy(lag2.mods[[1]])[2,],
                  tidy(lag2.mods[[3]])[2,],
                  tidy(lag2.mods[[7]])[2,],
                  tidy(lag2.mods[[8]])[2,],
                  tidy(lag2.mods[[9]])[2,]) %>%
  mutate(lag = 2,
         std.error = c(ses[[1]][2],
                       ses[[3]][2],
                       ses[[7]][2],
                       ses[[8]][2],
                       ses[[9]][2]))


long.lags.3 <- working.data.all.tsmos %>%
  group_by(ccode) %>%
  arrange(year) %>%
  mutate_at(vars(tsmo.wtd.polyarchy:tsmo.wtd.suffr), ~ lag(.,n = 2L)) %>%
  ungroup()

attach(long.lags.3)
lag3.mods <- tsmomodel(dvform = "lead",
                       ivform = "wtd",
                       controlvars = c("ioscore","dem.perf","cw","log_gdppc","neighbor","log.tsmo.mems"))
 ses <- lapply(lag3.mods, function(x) {
    sqrt(diag(cluster.vcov(x, ccode)))
  })
detach(long.lags.3)

lag3 <- bind_rows(tidy(lag3.mods[[1]])[2,],
                  tidy(lag3.mods[[3]])[2,],
                  tidy(lag3.mods[[7]])[2,],
                  tidy(lag3.mods[[8]])[2,],
                  tidy(lag3.mods[[9]])[2,]) %>%
  mutate(lag = 3,
         std.error = c(ses[[1]][2],
                       ses[[3]][2],
                       ses[[7]][2],
                       ses[[8]][2],
                       ses[[9]][2]))

long.lags.4 <- working.data.all.tsmos %>%
  group_by(ccode) %>%
  arrange(year) %>%
  mutate_at(vars(tsmo.wtd.polyarchy:tsmo.wtd.suffr), ~ lag(.,n = 3L)) %>%
  ungroup()

attach(long.lags.4)
lag4.mods <- tsmomodel(dvform = "lead",
                       ivform = "wtd",
                       controlvars = c("ioscore","dem.perf","cw","log_gdppc","neighbor","log.tsmo.mems"))
 ses <- lapply(lag4.mods, function(x) {
    sqrt(diag(cluster.vcov(x, ccode)))
  })
detach(long.lags.4)

lag4 <- bind_rows(tidy(lag4.mods[[1]])[2,],
                  tidy(lag4.mods[[3]])[2,],
                  tidy(lag4.mods[[7]])[2,],
                  tidy(lag4.mods[[8]])[2,],
                  tidy(lag4.mods[[9]])[2,]) %>%
  mutate(lag = 4,
         std.error = c(ses[[1]][2],
                       ses[[3]][2],
                       ses[[7]][2],
                       ses[[8]][2],
                       ses[[9]][2]))

long.lags.5 <- working.data.all.tsmos %>%
  group_by(ccode) %>%
  arrange(year) %>%
  mutate_at(vars(tsmo.wtd.polyarchy:tsmo.wtd.suffr), ~ lag(.,n = 4L)) %>%
  ungroup()

attach(long.lags.5)
lag5.mods <- tsmomodel(dvform = "lead",
                       ivform = "wtd",
                       controlvars = c("ioscore","dem.perf","cw","log_gdppc","neighbor","log.tsmo.mems"))
 ses <- lapply(lag5.mods, function(x) {
    sqrt(diag(cluster.vcov(x, ccode)))
  })
detach(long.lags.5)

lag5 <- bind_rows(tidy(lag5.mods[[1]])[2,],
                  tidy(lag5.mods[[3]])[2,],
                  tidy(lag5.mods[[7]])[2,],
                  tidy(lag5.mods[[8]])[2,],
                  tidy(lag5.mods[[9]])[2,]) %>%
  mutate(lag = 5,
         std.error = c(ses[[1]][2],
                       ses[[3]][2],
                       ses[[7]][2],
                       ses[[8]][2],
                       ses[[9]][2]))

long.lags.6 <- working.data.all.tsmos %>%
  group_by(ccode) %>%
  arrange(year) %>%
  mutate_at(vars(tsmo.wtd.polyarchy:tsmo.wtd.suffr), ~ lag(.,n = 5L)) %>%
  ungroup()

attach(long.lags.6)
lag6.mods <- tsmomodel(dvform = "lead",
                       ivform = "wtd",
                       controlvars = c("ioscore","dem.perf","cw","log_gdppc","neighbor","log.tsmo.mems"))
 ses <- lapply(lag6.mods, function(x) {
    sqrt(diag(cluster.vcov(x, ccode)))
  })
detach(long.lags.6)

lag6 <- bind_rows(tidy(lag6.mods[[1]])[2,],
                  tidy(lag6.mods[[3]])[2,],
                  tidy(lag6.mods[[7]])[2,],
                  tidy(lag6.mods[[8]])[2,],
                  tidy(lag6.mods[[9]])[2,]) %>%
  mutate(lag = 6,
         std.error = c(ses[[1]][2],
                       ses[[3]][2],
                       ses[[7]][2],
                       ses[[8]][2],
                       ses[[9]][2]))

long.lags.7 <- working.data.all.tsmos %>%
  group_by(ccode) %>%
  arrange(year) %>%
  mutate_at(vars(tsmo.wtd.polyarchy:tsmo.wtd.suffr), ~ lag(.,n = 6L)) %>%
  ungroup()

attach(long.lags.7)
lag7.mods <- tsmomodel(dvform = "lead",
                       ivform = "wtd",
                       controlvars = c("ioscore","dem.perf","cw","log_gdppc","neighbor","log.tsmo.mems"))
 ses <- lapply(lag7.mods, function(x) {
    sqrt(diag(cluster.vcov(x, ccode)))
  })
detach(long.lags.7)

lag7 <- bind_rows(tidy(lag7.mods[[1]])[2,],
                  tidy(lag7.mods[[3]])[2,],
                  tidy(lag7.mods[[7]])[2,],
                  tidy(lag7.mods[[8]])[2,],
                  tidy(lag7.mods[[9]])[2,]) %>%
  mutate(lag = 7,
         std.error = c(ses[[1]][2],
                       ses[[3]][2],
                       ses[[7]][2],
                       ses[[8]][2],
                       ses[[9]][2]))

long.lags.8 <- working.data.all.tsmos %>%
  group_by(ccode) %>%
  arrange(year) %>%
  mutate_at(vars(tsmo.wtd.polyarchy:tsmo.wtd.suffr), ~ lag(.,n = 7L)) %>%
  ungroup()

attach(long.lags.8)
lag8.mods <- tsmomodel(dvform = "wtd",
                       ivform = "lead",
                       controlvars = c("ioscore","dem.perf","cw","log_gdppc","neighbor","log.tsmo.mems"))
 ses <- lapply(lag8.mods, function(x) {
    sqrt(diag(cluster.vcov(x, ccode)))
  })
detach(long.lags.8)

lag8 <- bind_rows(tidy(lag8.mods[[1]])[2,],
                  tidy(lag8.mods[[3]])[2,],
                  tidy(lag8.mods[[7]])[2,],
                  tidy(lag8.mods[[8]])[2,],
                  tidy(lag8.mods[[9]])[2,]) %>%
  mutate(lag = 8,
         std.error = c(ses[[1]][2],
                       ses[[3]][2],
                       ses[[7]][2],
                       ses[[8]][2],
                       ses[[9]][2]))

long.lags.9 <- working.data.all.tsmos %>%
  group_by(ccode) %>%
  arrange(year) %>%
  mutate_at(vars(tsmo.wtd.polyarchy:tsmo.wtd.suffr), ~ lag(.,n = 8L)) %>%
  ungroup()

attach(long.lags.9)
lag9.mods <- tsmomodel(dvform = "wtd",
                       ivform = "lead",
                       controlvars = c("ioscore","dem.perf","cw","log_gdppc","neighbor","log.tsmo.mems"))
 ses <- lapply(lag9.mods, function(x) {
    sqrt(diag(cluster.vcov(x, ccode)))
  })
detach(long.lags.9)

lag9 <- bind_rows(tidy(lag9.mods[[1]])[2,],
                  tidy(lag9.mods[[3]])[2,],
                  tidy(lag9.mods[[7]])[2,],
                  tidy(lag9.mods[[8]])[2,],
                  tidy(lag9.mods[[9]])[2,]) %>%
  mutate(lag = 9,
         std.error = c(ses[[1]][2],
                       ses[[3]][2],
                       ses[[7]][2],
                       ses[[8]][2],
                       ses[[9]][2]))

long.lags.10 <- working.data.all.tsmos %>%
  group_by(ccode) %>%
  arrange(year) %>%
  mutate_at(vars(tsmo.wtd.polyarchy:tsmo.wtd.suffr), ~ lag(.,n = 10L)) %>%
  ungroup()

attach(long.lags.10)
lag10.mods <- tsmomodel(dvform = "wtd",
                       ivform = "lead",
                       controlvars = c("ioscore","dem.perf","cw","log_gdppc","neighbor","log.tsmo.mems"))
 ses <- lapply(lag10.mods, function(x) {
    sqrt(diag(cluster.vcov(x, ccode)))
  })
detach(long.lags.10)

lag10 <- bind_rows(tidy(lag10.mods[[1]])[2,],
                   tidy(lag10.mods[[3]])[2,],
                   tidy(lag10.mods[[7]])[2,],
                   tidy(lag10.mods[[8]])[2,],
                   tidy(lag10.mods[[9]])[2,]) %>%
  mutate(lag = 10,
         std.error = c(ses[[1]][2],
                       ses[[3]][2],
                       ses[[7]][2],
                       ses[[8]][2],
                       ses[[9]][2]))

all_lags <- bind_rows(lag1,lag2,lag3,lag4,lag5,lag6,lag7,lag8,lag9,lag10) %>%
  mutate(term = case_when(term == "tsmo.wtd.polyarchy" ~ "Polyarchy",
                          term == "tsmo.wtd.fr.assoc" ~ "Free Association",
                          term == "tsmo.wtd.frefair" ~ "Free Elections",
                          term == "tsmo.wtd.fre.exp" ~ "Free Expression",
                          TRUE ~ "Participatory Democracy"))

lag.plot <- ggplot(all_lags,aes(lag,estimate)) +
             geom_point() +
             geom_errorbar(aes(ymin=estimate - (1.96*std.error),ymax = estimate + (1.96*std.error)),width = 0.1) +
  geom_hline(yintercept = 0) +
  theme_minimal() +
  labs(x = "Years Lagged",y = "Coefficient Estimate",title = "TSMO Score Long-Term Robustness") +
  scale_x_continuous(breaks = c(2,4,6,8,10)) +
  facet_wrap(~ term,ncol = 2)
lag.plot

```
